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We investigate theoretically and experimentally the statistical properties of the inhomegeneous 
order-parameter distribution (OPD) at the verge of the superconducting-insulator transition (SIT). 
We find within two prototype fermionic and bosonic models for disordered superconductors that 
one can identify an universal rescaling of the OPD. By performing scanning-tunneling microscopy 
experiments in three samples of NbN with increasing disorder we show that such a rescaling describes 
also with an excellent accuracy the experimental data. These results can provide a breakthrough in 
our understanding of the SIT. 



PACS numbers: 74.20.Mn 71.30.+h 74.20.-z 71.55.Jv 

The interplay between disorder and superconductivity 
represents a typical example of emerging complex behav- 
ior in the presence of competing mechanisms. Indeed, 
while the former leads to localization of the electrons, 
leading to insulating-like transport, the latter favors the 
formation of a macroscopic coherent electronic state able 
to substain a dissipationless current. While at moder- 
ate disorder level the pairing mechanism persists almost 
unchanged [L, as disorder increases the superconducting 
(SC) critical temperature T c decreases and ultimately 
a full insulating state is reached. The most interesting 
case occurs when the superconductor-insulator transition 
(SIT) is somehow direct, i.e. without an intermediate 
bad-metallic state. Indeed, in this situation one can ex- 
pect a persistence of SC correlations in the insulating 
state and conversely a precursor effect of the insulating 
order on the SC side[5J[3]. 

In the last few years considerable theoretical and ex- 
perimental advances have been made to outline such a 
scenario on solid grounds. In particular, new insights 
have been offered by experiments of scanning tunnel- 
ing microscopy [4- -8J, that have access to the local den- 
sity of states (DOS) of homogeneously strongly dis- 
ordered superconductors. The most striking features 
are the emergence of an intrinsic mesoscopic inhomo- 
geneity in the local SC properties, and the occurrence 
of a large scale Ap ^> T c for the DOS suppression, 
that persists well above T c . These effects are under- 
stood qualitatively by using prototype models of disor- 
dered superconductors [SJ, that can be based either on a 
fermionic [TTjMTfj] or bosonic [TTHTU] description of the rel- 
evant degrees of freedom. In the former case it has been 
demonstrated the contribution of strong disorder to man- 
tain the scale Ap finite at the SIT, driven essentially by 
a loss of SC phase coherence. In the latter case the focus 
has been put on a SIT driven mainly via a localization of 



preformed pairs, due to quantum fluctuations associated 
to the random local energies. Within this scenario it has 
been also predicted|18l HTj] a non-self-averaging behav- 
ior for the inhomogeneous order-parameter distribution, 
with a characteristic power-law decay which implies a di- 
verging mean in the absence of the upper physical cutoff. 

Despite some valuable attempts [51 [TJ] to establish a 
link between theoretical predictions and experiments, it 
has not been established yet which characteristic signa- 
tures of the SIT allow one for a convincing quantita- 
tive comparison between theory and experiments. The 
present work aims at filling this gap, by establishing at 
the same time a bridge between the two lines of theoret- 
ical investigations mentioned above. More specifically, 
in analogy with Refs. [6j [HI [19], we shall focus on the 
properties of the order-parameter distribution (OPD), in 
order to demonstrate the emergence of universal scaling 
properties at the verge of the SIT. We show that both 
within fermionic and bosonic models the OPD evolves at 
strong disorder not only with a diverging mean of the 
order parameter (OP), but also with a diverging width 
of the logarithm of the OP. This suggests an universal 
scaling of the OPD in the SC phase that covers not 
only the theoretical results obtained with different ap- 
proaches, but more remarkably is in excellent agreement 
with experimental data in three different samples of dis- 
ordered NbN films. We notice that such OPD differs 
from the one suggested in Refs. (T5J ITS] , Their result 
was obtained by a mapping into the directed-polymer 
model on the Cayley tree, where the number of neigh- 
bors grows exponentially with the distance, in contrast 
to ordinary finite-dimensional lattices. In this respect 
our result establish instead the relevance of the directed- 
polymer physics in finite dimension [5T] [55J for the SIT. 

The first prototype fermionic model for a disordered 
superconductor that we will analyze is the Hubbard 
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Figure 1: (Color online) Disorder-dependence of the OPD within(a) C-CMF, (b) BdG and (c) 2D-CMF. The probability for 
each S scales as the color code shown on the right of each panel. The maximum of the OPD is located approximately at the 
typical OP, Styp, whose g dependence is shown with a continous line in the main panels. The insets show explicitly the In S 
dependence of the OPD for selected representative g values in the superconducting phase, marked by vertical bars in the main 
panels. Notice that in C-CMF a power-law behavior P(lnS) ~ S -0 ' 7 sets in for large OP values, 5t yp S g/K as predicted 
in [19] (see the brown dashed line in the inset). Instead within BdG and 2D-CMF one observes the formation of strongly 
asymmetric distributions with large tails extending towards small 5 values. The parameters are the following (see text): (a) 
L = 10, (b) \U\ = 5, (n) = 0.875, (c) L = 1000. 



model with random on-site energies: 

H = -t {c\ a c ]a +h.c)+y^\V i -^)n ia -\U\y^ j n^n i ±. 

(i) 

Here c\ a (c; CT ) is the creation (destruction) operator for 
an electron with spin a on a site of a square lattice 
with lattice spacing a = 1, t = 1 is the nearest-neighbor 
hopping, \U\ is the pairing interaction, n.j CT = cj^Cjo-, 
and /i is the chemical potential. The on-site potentials 
Vi are independent quenched random variables box dis- 
tributed between — Vo and Vo, with Vo denoting the dis- 
order amplitude. We will investigate the model ([l} by 
means of Bogoliubov-de Gennes (BdG) mean field the- 
ory pTj [TBI E2], by allowing for spatial fluctuations of 
the SC OP Aj = \U\(ci^Ci^) . Even though within such 
a mean-field approach one cannot describe the SIT, it 
captures already several features of strongly-disordered 
superconductors [TTJ [T31 [21], as the emergence of spatial 
inhomogeneity of the OP and of a large DOS suppres- 
sion due to the interplay between superconductivity and 
disorder. More recently it has been shownQJj] also that 
at strong disorder the SC current follows a non-trivial 
percolative pattern, reminiscent of the glassy behavior 
suggested by the analysis of Ref . [HI [19] . Here we want 
to explore more deeply such analogy, by comparing the 
results obtained within the model ([IJ with the effective 
bosonic model investigated in Refs. [TBI [El], i- e - the Ising 
spin model in a random transverse field: 

i <ij> 

The model (12]) puts the main emphasis on the compe- 



tition between local pairing and single-particle localiza- 
tion that is realized at strong disorder [10 . Thus, in this 
language a state with of = ±1 corresponds to a site 
occupied or unoccupied by a Cooper pair, while the su- 
perconducting phase corresponds to the existence of a 
spontaneous magnetization in the x direction. The on- 
site energies which play the role of the transverse field 
in the spin language, are independent quenched random 
variables which are chosen here to be box distributed be- 
tween — 1 and 1. The insulating phase corresponds thus 
to spins aligned along the z axis. The connection be- 
tween the two approaches is suggested by a mapping be- 
tween a disordered superconductor and a random (XY) 
spin mo del [TO], where the low-energy degrees of free- 
dom are pairs that can hop from one site to another. 
In the strong-coupling regime U ^> t one can also show 
[2"5] that in the clean case the effective spin coupling is 
g = t 2 /VoU . Notice that the model ^ lacks completely 
the XY symmetry that allows for the existence of the 
Goldstone mode, i.e. for the SC phase fluctuations. In 
this respect a comparison with the BdG approach, which 
also neglects phase fluctuations, is meaningful. In the 
following we will show that this approximation is enough 
to describe the anomalous effects of the OPD distribu- 
tion at the verge of the SIT, since this physics is driven 
mainly by the competition between pair hopping and site 
localization. This does not exclude of course that at the 
SIT phase fluctuations will lead to additional remarkable 
effects, as discussed in Refs. [T3HT6] and suggested exper- 
imentally by measurements of the penetration depthjS]. 

The Hubbard model ([!]) has been investigated in a wide 
range of parameters and lattices of linear dimensions up 
to L = 25, with a large number of disorder configurations 
(up to 1920). As far as the Ising model ([2| is concerned, 



in addition to a standard mean-field approach (2D-MF) 
we used also a cavity- mean-field approximation [18H2Q], 
that allows one to include quantum effects which lead to 
a SIT for a finite value of the coupling g. To emphasize 
also the role of the underlying lattice we will show results 
for two different lattice structures: the Caley tree with 
K = 3 (C-CMF), discussed in Refs. OH EH], and the 2D 
lattice (2D-CMF) discussed in Ref. [21], where the cavity 
approximation is implemented via a non-linear transfer 
approach. For the random Ising model ([2| defined on 
a Caylee tree of connectivity K, the C-CMF approxi- 
mation [TBI EH] describes the spin j by the local Hamil- 
tonian ff? MF = - a^^tiK) where (a%) is 

the magnetization at site k due to the rest of the spins 
in absence of j and g = Kg is the coupling parameter 
considered in this context and which we will continue to 
denote g in the following. By defining the cavity mean 
field Bj = 4= SfcLi(°1c)> one obtains at zero temperature 
the CMF recursion relation: 

K B 

The 2D-CMF approach propagates instead equation |3]) 
with K — 2 along the diagonals of a 2D square lattice 
(see [2TJ for more details). 

In Fig. [T]we show the evolution with the coupling pa- 
rameter g of the OPD P(lnS) for C-CMF, BdG, and 
2D-CMF. Here Si — Bi/g is the normalized cavity field, 
which corresponds within BdG to a coarse-graining of the 
A, over nearest neighbors, Si = \ J2k=i TT7T • ^ ne P ror> 
ability at each In S value for the given disorder strength 
is represented in a color plot, where the maximum of the 
distribution is located approximately at the typical value 
of the OP, S typ = exp(Ih~5). In 2D-CMF and C-CMF 
the SIT transition occurs at g c « 0.22 and g c « 0.11, 
respectively, while in BdG the system remains supercon- 
ducting right up to g c = 0. As one can see in the insets 
of Fig.EjJi, where -P(ln£>) is reported for some representa- 
tive g values, for C-CMF we recover the expected power- 
law decay P(lnS) ~ S~ m with the universal (disorder- 
independent) exponent m = 1 — eg c ~ 0.7 predicted in 
Refs. CES] [JJ5]. However, such a power-law behavior is 
absent in the BdG and 2D-CMF results, where instead 
P(\nS) appears to be dominated by the low-field val- 
ues, and to be strongly disorder-dependent. We notice 
that such a discrepancy between C-CMF and 2D results 
cannot be attributed to the method itself: indeed, as we 
discuss in the Supplementary Material, in ID CMF and 
BdG give both p(S) ~ S a with a — > — 1 + when g — > g c 
which is in agreement with the exact critical behavior of 
the Ising model (Eg) [20"ll22]. 

Such a distinction between C-CMF from one side and 
2D results from the other can be made more quantitative 
by introducing as a scaling variable the logarithm of the 
OP, normalized to its variance a s = In 2 S — In S . In- 




Fi gure 2: (Color online) (a) Evolution of the typical OP Styp 
and of the distribution width as with increasing disorder (i.e. 
decreasing S tyP ) for 2D-CMF, 2D-MF, BdG and C-CMF. No- 
tice that while within C-CMF as saturates for increasing dis- 
order the 2D results show all an increasing as ■ We also report 
the three points corresponding to the NbN samples analized 
in Fig.[3]below. (b) Rescaling of the OPD with respect to St yp 
and as- All the 2D results collapse in a single curve, well fitted 
by the Tracy- Widom distribution, while the C-CMF results 
follow a different behavior. The parameters are the following: 
BdG (a), \U\ = 9, (n) = 0.3, L = 25, g = 0.1; BdG (b), 
\U\ = 5, (n) = 0.875, L = 25, g = 0.08; 2D-CMF, L = 1000, 
g = 0.4; 2D-MF, L = 120, g = 0.2; C-CMF, L = 15, K = 3 
and g = 0.2. 



deed, as one can see in Fig. [2^,, when disorder increases 
S typ and ag scale in the same way in the 2D case, while 
within C-CMF as tends to saturate at strong disorder. 
This result hints to a remarkable property of the OPD, 
that becomes evident when the above data are rescaled 
with the variable IZs = (hi S — ln<Styp)/o\s- As shown in 
Fig. [2}d, provided that the coupling is small enough but 
still in the SC phase, all the data (except the ones for the 
C-CMF) collapse into a single curve, well approximated 
by the Tracy Widom distribution[26 , whose relevance in 
the insulating phase has been recently discussed in Refs. 
[2TJ [22] . Notice that such scaling holds for a wide range 
of parameters both within fermionic and bosonic models, 
showing the robustness of the present result. 

To compare these results with experiments, we per- 
formed scanning tunneling spectroscopy (STS) measure- 
ments at 500 mK on three NbN films with different lev- 
els of disorder, having T c ^1.68 K, 2.9 K and 6.1 K re- 
spectively. NbN seems to be the ideal candidate to in- 
vestigate the SIT, since it has been shown earlier that 
its T c decreases monotonically with increasing disorder, 
eventually giving rise to a non-SC state characterized by 
strong SC correlations [7]. The thickness of all the films 
was ^50 nm, which is much larger than the dirty-limit 
coherence length [25]. Details of sample deposition and 
characterization have been reported in Refs. [71 [311] . 
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Figure 3: (Color online) (a)-(b) Representative background- 
corrected tunneling spectra at 500 mK for two different lo- 
cations on a NbN film with T c ~1.65 K; the horizontal red 
line shows the normal-state conductance. The insets show the 
corresponding tunneling spectra before background subtrac- 
tion (blue line) along with the spatially-averaged spectrum at 
8K used to subtract the background (red line). The average 
peak height h = (hi + h,2)/2 measures the local order param- 
eter, whose spatial variation measured at 500 mK is shown 
in panels (c) (T c ~1.65 K ) and (d) (T c ~2.9 K) over a 200 
nm X 200 nm area. Panel (e) shows the OPD in the three 
samples in linear scale. The same data are reported in panel 
(f) using the rescaling adopted in Fig. |2jj. The theoretical 
curve corresponds to the Tracy-Widom distribution. 



For STS measurements, the samples were grown in-situ 
in a sputtering chamber connected to a home-built scan- 
ning tunneling microscope (STM)[57]. For each film, 
tunneling conductance (dl/dV vs V) was measured at 
each point on a 32 X 32 grid over an area of 200 X 
200 nm. All the spectra show a prominent dip in the 
conductance associated with the superconducting energy 
gap which rides over a broad V-shaped background aris- 
ing from Altshuler-Aronov type electron-electron inter- 
actions which extends up to high bias (inset Fig. 3a, 3b) . 
To isolate the feature associated with superconductivity 
from the background we follow the procedure adopted in 
Refs. [5J [7] : The spatially averaged tunneling spectrum 
recorded at 8 K where all superconducting correlations 
are destroyed is used to subtract the background from the 
spectra recorded at low temperatures. The background- 
corrected tunneling spectra do not show a significant vari- 
ation in the magnitude of the superconducting energy 
gaps|31). but they show a large distribution in the height 



of the coherence peaks: At some positions we observe 
a prominent coherence peak (Fig. 3a) whereas at other 
places the coherence peak is completely suppressed (Fig. 
3b). We take the average of the heights of the coherence 
peaks at positive and negative bias from the normal state 
conductance values h = ((hi + /i2)/2) as the measure of 
the local OP Aj for our system [HI US]- Figures 3c and 
3d show the spatial variation of the OP in the form of 
an intensity plot for the samples with T c ~ 1.65 K and 
T c ~ 2.9 K, respectively. We observe smooth variation 
of Aj over length scales of few tens of nanometers, lead- 
ing to the OPD for the three samples reported in Fig. [3j3 
in linear scale. Here we defined S = h/M&x[h] to keep 
the same overall normalization used for the theoretical 
results. As disorder increases one observes a steady de- 
crease of the OPD maximum along with a widening of 
the OPD, similar to the one reported in TiN samples in 
Ref. [6]. This can be further quantified by computing 
S^yp and the variance of the log a e g P , which are found to 
follow the same trend as the theoretical 2D results (see 
Fig. [2^,) . However the most striking result is that by in- 
troducing the scaling variable TZs = (In S — lniStyp)/cs 
defined above the three experimental OPD distributions 
collapse into a single universal curve, despite their appar- 
ent difference when plotted in linear scale. In addition, 
the agreement with the universal Tracy-Widom distribu- 
tion found in finite dimensions is also very good, even 
though it could be further improved at positive 1Z values 
to fully capture the striking universality of the experi- 
mental data. 



In summary we have shown both theoretically and ex- 
perimentally that the SC state at the verge of the SIT 
transition is characterized by an universal behavior of the 
OPD. The relevant scaling variable is the logarithm of the 
OP normalized to its variance. The latter diverges as the 
SIT is approached, unless the problem is studied on an 
infinite-dimensional lattice as the Caley tree, explaining 
the lack of such universality within the C-CMF [T51 fT§] . 
The universal OPD shares a pronounced similarity with 
the Tracy-Widom distribution, whose role in the disor- 
dered phase of the random Ising model has been recently 
discussed within the mapping into the directed-polymer 
model in finite dimensions [3TJ [52]. Within such a map- 
ping additional predictions have been made, as e.g. the 
divergence of the dynamical critical exponent as the SIT 
is approached. This could be tested experimentally by 
the critical scaling of the superconducting fluctuations 
at T c , as done recently in other systems [35]. While the 
critical properties of real systems at the SIT should ulti- 
mately belong to the XY universality class, at interme- 
diate disorder further experimental and theoretical in- 
vestigation of these predictions will further clarify the 
relevance of the directed-polymer physics on the SIT. 
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